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Abstract. Several stochastic simulation algorithms (SSAs) have been recently 
proposed for modelling reaction-diffusion processes in cellular and molecular biology. In 
this paper, two commonly used SSAs are studied. The first SSA is an on-lattice model 
described by the reaction-diffusion master equation. The second SSA is an off-lattice 
model based on the simulation of Brownian motion of individual molecules and their 
reactive collisions. In both cases, it is shown that the commonly used implementation 
of bimolecular reactions (i.e. the reactions of the form A+B — » C, or A+A — > C) might 
lead to incorrect results. Improvements of both SSAs are suggested which overcome the 
difficulties highlighted. In particular, a formula is presented for the smallest possible 
compartment size (lattice spacing) which can be correctly implemented in the first 
model. This implementation uses a new formula for the rate of bimolecular reactions 
per compartment (lattice site). 
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1. Introduction 



Many cellular and subcellular biological processes can be described in terms of diffusing 
and chemically reacting species (e.g. enzymes) [TJ [11]. A traditional approach to 
the mathematical modelling of such reaction-diffusion processes is to describe each 
(bio)chemical species by its (spatially-dependent) concentration. The time evolution 
of concentrations is then modelled by a system of partial differential equations (PDEs) 
[24]. Many mathematical and computational methods have been developed over the 
last century for solving and analyzing PDEs [291 S3] , which makes PDE-based modelling 
attractive. However, it has serious limitations when applied to biological systems. There 
may be relatively few numbers of some chemical species; for example, often only one 
or two niRNA molecules of a particular gene are present in the cell pp. In such cases 
we cannot even properly define spatially-dependent concentration profiles^ and PDE- 

* The macroscopic concentration of molecules at a given point in the space is defined as the number of 
molecules in a neighbourhood of this point divided by the volume of the neighbourhood. In particular, 
the neighbourhood must be chosen large enough to contain a lot of molecules. This is clearly not 
possible if there are only few molecules present in the system. 
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based models cannot be used. The appropriate quantities to describe the system are not 
concentrations, but numbers and positions of molecules of the chemical species involved. 

In recent years, several stochastic simulation algorithms (SSAs) have been proposed 
to model the time evolution of molecular numbers [20l [21 [TTJ. They provide a more 
detailed and precise picture than deterministic PDE-based models. They typically 
give the same results for simple systems involving zero-order and first-order chemical 
reactions (for example, linear degradation or conversion). However, the situation is 
more delicate whenever some chemical species are subject to bimolecular (second-order) 
reactions, or the system under study includes reactive boundaries (for example, a cellular 
membrane with receptors). Reactive boundaries were studied in our previous paper 1 1 1 j . 
where we systematically investigated four different SSAs for reaction-diffusion processes 
which had been proposed in the literature. We showed that one would obtain incorrect 
results if the computer implementation of reactive boundaries is not handled with care. 
In particular, what seems on the face of it the same boundary condition leads to different 
results when applied to different SSAs. To fix this problem, we derived formulae giving 
the correct relation between experimentally measurable characteristics and parameters 
of the computer implementation of boundary conditions for all four SSAs [11] . A 
generalization of one of these formulae to anisotropic diffusion tensors was recently 
given in [26J. 

In this paper, we focus on modelling bimolecular reactions, i.e. chemical reactions of 
the form A+B — > C or A+A — > C. We investigate two commonly used reaction-diffusion 
SSAs which have been previously implemented in reaction- diffusion software packages 
MesoRD [20] and Smoldyn [2]. The first reaction-diffusion SSA is based on dividing the 
computational domain into artificially well-mixed compartments and postulating that 
only molecules which are within the same compartment can react. Diffusion is then 
modelled as jumps between the neighbouring compartments. This approach can be 
mathematically described by the reaction-diffusion master equation [221 El EE] and was 
recently implemented in the mesoscopic reaction-diffusion simulator MesoRD [20]. In 
order to use this method, we have to choose an appropriate compartment size. On one 
hand, the compartment size must be chosen small enough so that the spatial variation 
in the concentration profiles can be captured with a desired resolution. The situation 
is analogous to solving PDEs numerically by a finite difference method. In order to 
solve PDEs with the desired accuracy, we need to choose a sufficiently fine mesh for 
discretization. On the other hand, we will see in Section 13.11 that the compartment size 
cannot be chosen arbitrarily small. The analogy with PDEs fails here. Unlike in the 
case of PDEs (for which we get a more accurate solution by using a finer discretization), 
there is a limit on the compartment size from below. In Section 13.11 we will show that 
the error of the computation increases as the compartment size decreases. In Section 
14.11 we present the formula for the smallest compartment size (which can be simulated 
by this approach) and propose an improved SSA which minimizes the simulation error, 
by modifying the reaction rate per compartment. 

The second SSA studied in this paper is based on Brownian motion of individual 
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molecules. In its classical formulation [27J , it is postulated that two molecules (which are 
subject to a bimolecular reaction) react whenever they are within a specified distance 
(reaction radius) from each other. A variant of this method was recently implemented in 
the software package Smoldyn [2] . One disadvantage of this approach is that the reaction 
radius is, for typical values of the bimolecular rate constant and diffusion coefficient, 
unrealistically small compared to the size of individual molecules. In Section 14.21 we 
propose an improved SSA to overcome this difficulty. It is based on the assumption 
that two molecules react with the rate A whenever they are within the distance ~q. The 
formula relating A, ~q and the simulation time step with the experimentally measurable 
reaction rate constant is derived. This formula is used for developing a more realistic 
SSA for reaction-diffusion processes. 

The paper is organized as follows. In Section [2], we present illustrative examples 
which will be used to demonstrate the results of the paper. In Section [3j we present 
both react ion- diffusion SSAs and summarize their major disadvantages. In Section HJ 
we present modified algorithms which are able to overcome the problems highlighted 
in Section [31 To make this paper accessible to non-mathematicians, Section [3] only 
contains the description of improved algorithms and formulae, together with the results 
of illustrative computations. The mathematical derivation of the formulae presented 
and the justification of the modified algorithms are given in Appendices. We finish with 
a discussion and conclusions in Section [51 

2. Bimolecular reactions - two model problems 

A bimolecular reaction is a chemical reaction involving two reacting molecules. 
Examples include 

A + B C + D, A + A C, or A + B B, 

where the capital letters A, B, C stand for chemical species and k is the reaction 
rate constant, expressed in units of volume over time. From the modelling point 
of view, it is useful to divide bimolecular reactions into two classes, heteroreactions 
and homoreactions. The term heteroreaction will be used for the bimolecular reaction 
between molecules of two different chemical species (for example, heterodimerization 
A + B —>■ C or catalytic degradation A + B — > B) . The bimolecular reaction between 
two molecules of the same chemical species (for example, homodimerization A + A — > C) 
will be called the homoreaction in what follows. In this section, we introduce two simple 
chemical systems which will be used to illustrate the results in the paper. The first 
model will include a heteroreaction (catalytic degradation) and the second model a 
homoreaction (homodimerization). More complicated examples are discussed later in 
Section [5j 
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2.1. A heteroreaction example 

Let us consider chemical species A and B in a container of volume v which are subject 
to the following two chemical reactions 

A + B 5, A (1) 

The first reaction is the degradation of A catalyzed by B. We couple it with the second 
reaction which represents the production of molecules of A with the ratjf] k 2 v. Since 
the number of B molecules is preserved in the chemical reactions (pQ), the dynamics of 
the model (CD) is simple: some molecules of A are produced by the second reaction and 
some are destroyed by the first reaction. Thus, after an initial transient behaviour, the 
number of A molecules fluctuates around its equilibrium value. 

Let us consider first the case when the chemical system ([1]) is well-stirred. Then the 
probability of an ocurrence of the bimolecular reaction is proportial to the number of 
available pairs of reactants. Let us define the propensity functions of chemical reactions 
©by 

a x {t) = A(t)B(t)h/u, a 2 (t) = k 2 u (2) 

where A(t) is the number of molecules of A at time t, B(t) is the number of molecules 
of B at time t, and v is the system volume. Then the probability that the z-th reaction 
occurs in the infinitesimally small time interval [t,t + dt) is equal to cti(t) dt, i = 1,2. 
Note that any heteroreaction A + B — > ■ ■ ■ has a propensity function equal to cti(t), 
while the propensity function of homoreactions differs; this is the reason why we discuss 
them separately. 

The chance of occurrence of each reaction is given by the corresponding propensity 
function (TSJ). If the first chemical reaction occurs, then one molecule of A is removed 
from the system; if the second chemical reaction takes places, then one molecule of A 
is added to the system. Given the values of the rate constants and the initial numbers 
of molecules of A and B, the stochastic model of ([1]) is uniquely specified and can be 
simulated by the Gillespie SSA [TFJ [18]. In Figure [H(a), we plot A(t), computed by the 
Gillespie SSA, as the solid line, using parameter values k\/u = 0.2 sec -1 , k 2 v = 1 sec -1 , 
A(0) = 5 and -B(O) = 1. As expected the number of molecules of A fluctuates around 
the average value which is 5 for our parameters. The nature of these fluctuations can 
be summarized in terms of the stationary distribution. To compute it, we record the 
values of A(t) at equal time intervals and create a histogram of the recorded values. 
Dividing the histogram by the total number of recordings, we obtain distribution 4>(n) 
which is plotted in Figure Wljo) as the grey histogram. Thus, <fi{n) is the probability 
that there are n molecules of A in the system, provided that the system is observed 
for a sufficiently long time. Note that since the Gillespie SSA makes use of random 
numbers to compute the time evolution of the system, the computed A(t) depends on 
a particular realization of the algorithm. Repeating the computation (with a different 

* Let us note that ki (the rate constant of the zero-order reaction) has physical dimension of units per 
volume per time. Consequently, k,2V is expressed in units per time. 
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Figure 1. Stochastic simulation of the system of chemical reactions (TI|) for k\/v = 
0.2 sec -1 , kiv — lsec -1 and one molecule of B in the system, (a) A(t) given by one 
realization of the Gillespie SSA (solid line) for the initial value A(0) = 5. The average 
value of A is plotted as the dashed line, (b) Stationary distribution 4>{n) obtained by 
long time simulation of the Gillespie SSA (grey histogram) and by formula (TjJ (circles). 



set of random numbers), we will obtain a different time evolution than the one plotted 
in Figure [TJa). However, the stationary distribution <p{n) is uniquely determined by the 
values of the rate constants ki, k 2 and the number of molecules of B in the container. 
It can be shown (see Appendix A ) that 0(n) is the Poisson distribution 

k 2 u^ 



(j){n) = — 



77! 



h B 



cxp 



kiB 



(3) 



where B is the (constant) number of molecules of B in the container. The results 
given by the formula ([3]) are plotted in Figure QJb) as the circles. We confirm that 
the stationary distribution 4>(n) is indeed given by (jBJ)- In what follows, we will use the 
stationary distributions of model problems to study the limitations of different stochastic 
react ion- diffusion methods. The average number of molecules of A in the container is 
given by (see Appendix A ) 



M, 



k 2 v 1 



(4) 



Using the parameter values of Figure [TJ we obtain M s = 5. This number is plotted in 
Figure Ufa) at the dashed line. The variance of the Poisson distribution ([3]) is equal to 
its mean M„. 



2.2. A homoreaction example 

Let us consider chemical species A and B in a container of volume v which are subject 
to the following two chemical reactions 

fci i-> a 



A + A 



B. 



A. 



(5) 
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The first reaction describes the dimerization of the chemical A with the rate constant k\. 
We couple it with the second reaction which represents the production of the chemical 
A with the rate constant k 2 . This reaction has been already studied in the example 
(PQ). In what follows, we will only be interested in the time evolution and stationary 
distribution of A. The dynamics of the model ([5]) is similar to ([I]): some molecules 
of A are produced by the second reaction and some are removed by the first reaction. 
Thus, Ait) fluctuates around its equilibrium value in a similar way as the trajectory in 
Figure Q] which has been computed for the heteroreaction example (CQ). If the reactor is 
well-stirred, the propensity functions of chemical reactions ([5]) are given by 

ai(t) = A(t)(A(t) - l)h/u, a 2 (t) = k 2 u, (6) 

i.e. the probability that the i-th reaction in (jSJ) occurs in the infinitesimally small 
time interval [t, t + dt) is cx.i{t) dt, i — 1, 2. If the first chemical reaction occurs, then two 
molecules of A are removed from the system; if the second chemical reaction takes place, 
then one molecule of A is added to the system. This uniquely specifies the stochastic 
model as a Markov chain which can be simulated by the Gillespie SSA. The stationary 



distribution of (jSJ) is given by (see Appendix A) 




0(w) = -T -M /^-i 2W-f— , n = 0,l,2,3,..., (7) 



where I n is the modified Bessel function of the first kind (see Glossary) and C is a 
positive constant given by the normalization ~Y^ n 4>{ n ) = ^he avera g e number of 



molecules of A in the container is given by (see Appendix A) 
M, = 





-i 

(8) 



In Appendix A , we show that the stationary number of molecules of A obtained 



by the standard ordinary differential equation (ODE) model of the chemical system 
(J5j is A s — v^Jk 2 j1h\. Note that this is not in general equal to M s given by 
formula (jSj). For example, in the following section, we use the parameter values 
h\jv = 0.2 sec -1 and k 2 v = 10 sec -1 . Then A s = 5 and M s = 5.13, i.e. the deterministic 
ODE does not provide the exact description of the stochastic mean. On the other 
hand, the difference between M s and A s is only 2.5% so that the ODE model gives 
a reasonable approximation of M s . However, the comparison of deterministic and 
stochastic modelling is not the focus of this paper. Our main goal is to highlight some 
limitations of current reaction-diffusion SSAs and present improvements of these models. 
Examples of chemical systems where the differences between the results of stochastic 
simulation and the corresponding deterministic approximation (ODEs) are significant 
can be found in [231 H3 El HI- 
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Figure 2. Domain [0, L] x [0, L] x [0, L] is divided into K 3 compartments of the 
volume h 3 — (L/K) 3 . The division of the domain for K = 8 is shown on the picture. 

3. Disadvantages of current SSAs for reaction-diffusion modelling 

In Section [21 we considered the illustrative examples (PQ) and (jSJ) as well-stirred chemical 
systems. In particular, their mathematical models did not explicitly involve a description 
of molecular diffusion. In this section, we couple chemical systems flTJ and (jSJ) with 
models of molecular diffusion and show key limitations of reaction-diffusion SSAs in 
the literature. In what follows, we assume that chemical species A and B diffuse with 
diffusion constants Da and Db, respectively, in the cubic container [0, L] x [0, L\ x [0, L\. 
We consider zero-flux (reflective) boundary conditions, i.e. whenever a molecule hits 
the boundary, it is reflected back. The implementation of more complicated (reactive) 
boundary conditions was studied in our previous paper [TT] . 

3.1. Compartment-based model 

In the compartment-based model, we divide the computational domain into small 
compartments which are assumed to be well-mixed. We postulate that only molecules 
in the same compartment can react according to bimolecular reactions. Diffusion is 
modelled as jumps of molecules between neighbouring compartments (201 122] . 

Let us consider the heteroreaction example (CQ). We divide the cubic domain 
[0, L] x [0, L\ x [0, L] into K 3 cubic compartments of volume h 3 where K > 1 and 
h = L/ K (see FigureEj). In general, the compartment-based model can be formulated for 
compartments which are not cubic and which are not of the same size (221 El] ■ However, 
for the purposes of this paper, it is sufficient to work with cubic compartments of the 
same size. They are the most natural choice, and are easy to implement computationally. 
Moreover, if the modeller does not use the uniform cubic mesh, it might be sometimes 
difficult to distinguish which results show a genuine property of the system and which are 
a consequence of the non-uniform mesh. Note that the uniform cubic mesh introduces 
an artificial anisotropy in the domain (e.g. compartments have different lengths along 
the side and along the diagonal). However, we will not explore potential consequences 



Stochastic reaction- diffusion processes 



8 



of this anisotropy in this paper. We will focus on the more fundamental problem: the 
appropriate choice of compartment size h. 

To precisely formulate the compartment-based SSA for the illustrative chemical 
system ([!]) in the reactor [0, L] x [0, L] x [0, L], we denote the compartments by indices 
from the set 

I ail — {(hj, k) | k are integers such that 1 < k < K}. 

Let Aijfc(i) (resp. Bij k (t)) be the number of molecules of the chemical species A (resp. 
B) in the (i,j,k)-th compartment at time t where (i,j,k) G I a ii- Diffusion is modelled 
as a jump process between neighbouring compartments. Let us define the set of possible 
directions of jumps 

E = {[1,0,0], [-1,0,0], [0,1,0], [0,-1,0], [0,0,1], [0,0,-1]}. 

For every k) G I a iu we also define 

Eijfc = {e G E | ((i, j, k) + e) G I a ii} , 

i.e Ejjfc is the set of possible directions of jumps from the (i,j, k)-th compartment. The 
compartment-based reaction-diffusion model can be written using the chemical reactions 
formalism as follows. We study a system of 2K 3 "chemical species" and Bij k for 
k) G l a u which are subject to the chemical reactions: 

A ijk + B ijk B ijk , A ijk , for (i, j, k) G I aU , (9) 

A ijk D -^> A ijk+e , for (i,j,k) G I a ii, e G E ijk , (10) 

B ijk B ijk+e , for (i,j,k) G hiu e G E iifc . (11) 

The chemical reactions (Q correspond to the chemical system ([1]) considered in each 
compartment. It is assumed that each compartment is effectively well-stirred. A 
molecules of A and a molecule of B which are in the same compartment can react 
according to the bimolecular reaction A + B — > B. On the other hand, two molecules 
in different compartments cannot react with each other. The propensity functions of 
reactions ([9]) are 

a ijki i(t) = A ijk (t)B ijk (t) h/h 3 , a ijk , 2 (t) = k 2 h 3 (12) 

where h 3 is the volume of the compartment. The propensity functions ( TT2l) can be 
derived using the same argument as replacing the volume v = L 3 of the whole reactor 
by the compartment volume h 3 . The reactions (Il0l - (fTTT) correspond to diffusive jumps 
between neighbouring compartments. The propensity functions of these "reactions" are 
equal to A ijk (t) D A /h 2 and B ijk {t) D B /h 2 . There are 2K 3 reactions in (jSJ), 6K 3 - 6K 2 
diffusion "reactions" in fflUj) and 6K 3 — 6K 2 diffusion "reactions" in ffTTl) because there 
are 6 possible directions to jump from each inner compartment and some directions are 
missing for boundary compartments. Thus we are able to formulate the compartment- 
based reaction-diffusion model as the chemical system of 2K 3 chemical species A^ k 
and Bij k which are subject to 1AK 3 — 12K 2 reactions (l9l)-(fTTT). The time evolution 
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Figure 3. (a) Heter vr •taction example ([7]). Stationary distribution ^>k(ti), defined by 
([H| and computed for K — 1,2, 20, and 100 by long time simulations of the Gillespie 
SSA. We use k\ = 0.2 /im 3 sec -1 , = 1 fixa sec" 1 , Da = Db = 1 (an 2 sec -1 , 
L = 1 ptm and Bq = 1. (b) Homoreaction example ([5j). Stationary distribution 4>k{it-) 
for K = 1,2, 20, and 200 computed by long time simulations of the Gillespie SSA. We 
use ki — 0.2 /im 3 sec -1 , &2 = 10 /im -3 sec -1 , Da = 1 /mi 2 sec -1 and L — 1 /im. 



of the chemical system (|9|)- (|TT|) can be simulated by the Gillespie SSA. It can be also 
equivalently described in terms of the reaction-diffusion master equation, which is given 
Appendix B| as equation (IB. 21) . The number of molecules of A in the whole container 



in 



[0, L) x [0, L) x [0, L) is given by 

A(t) = A ^(t)- 

(i,j,k)el a ii 

Let p n (t) be the probability that A{t) = n. Let 4>k (n) be the stationary distribution 
defined by (compare with (1A.2j) ) 

<p K {n) = lim p n (t). (13) 

t^oo 

Thus, <pK{ n ) is the probability that there are n molecules of A in the system, provided 
that the system is observed for long time. In particular, 0i(n) is equal to the stationary 
distribution (p{n) given by ([3]). Since production of A is homogeneous throughout the 
container, we would expect the distribution of A to be uniform in space, so that we 
should find that 4>k is independent of K. In Figure Eta), we present the stationary 
distributions (pxin) for K = 1,2,20 and 100 for the parameter values L = 1 /im, 
Da = D B = 1 /im 2 sec -1 , k\ = 0.2 /im 3 sec -1 , k 2 = 1 /im" 3 sec -1 and B = 1. 
In particular, we have the same rates for K — 1 as were used in Figure (TJ namely 
ki/L 3 = 0.2 sec -1 and k 2 L 3 = 1 sec -1 . Thus the stationary distribution <fii(n), plotted 
in Figure [3](a), is equal to the distribution <p{n) plotted in Figure [U^b). 

Increasing K (i.e. decreasing h), the stationary distribution 0^(n) moves to the 
right. The shift to the right is in agreement with the result of Isaacson [21] who showed 
that, in the theoretical limit h —>■ 0, the bimolecular reaction A + B — > is lost and 
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the compartment-based modelling of this reaction only recovers the diffusion process. 
In our case, we coupled the bimolecular reaction with the production of A molecules. 
The production rate per the whole domain is equal to 

a m,2{t) = k ^ = R3k 2h 3 = k 2 , 

(i,j,h)el a ii {i,j,k)el a u 

i.e. it is independent of fa. Thus, for small fa, the slower removal of A by the bimolecular 
reaction and unchanged production rate of A result in the shift of the stationary 
distribution 4>k(j>) to the right in Figure [3]^a), as K is increased (i.e. as fa = L/K 
is decreased). In Figure E](b), we present the results of a similar computation for the 
homoreaction example ([5]). In this case, the homodimerization A + A — > B is replaced 
by K 3 reactions A^ + A^ — ► for k) G I a ii- The propensity functions of these 
reactions are given by 

atijk,i(t) = A ijk [t)[A ijk [t) - 1) h/h 3 . (14) 

The production reaction and diffusion are treated as in (l9i)-( fT0j) . In Figure G^b), we 
present the stationary distributions 4>K[ n ) for four values of K. Notice that </>i(n) is 
equal to the stationary distribution (f>(n) given by (JTj). We observe the same phenomenon 
(shift of the histogram to the right) as in the case of the heteroreaction example. 

Although it is generally agreed in the literature that there is a bound on h from 
below [22j EI]; this bound is usually stated in the form fa k\j [Da + Db) or h ^> g 
where g is the binding radius for the molecular based Smoluchowski model - see equation 
(12~TT) and the discussion in Section 13.21 To satisfy these conditions in our particular 
example, we could simply choose h = L. However, the real importance of stochastic 
react ion- diffusion modelling is not in modelling of spatially homogeneous systems. If 
the system has some spatial variations (i.e. some parts of the computational domain are 
more preferred by molecules than the others) , then we obviously want to choose h small 
enough to capture the desired spatial resolution. This leads to the restriction on h from 
above, namely L ^> h. Thus it is suggested to choose h small (to satisfy L h) but 
not too small (to satisfy h ^> k\/{DA + Db)) [22J which leads to the important question 
what the optimal choice of h should be to get the most accurate results. In this paper, 
we propose a different route to this problem. In Section 14.11 we show that there exists 
a critical value h cr i t such that the propensity function of the compartment-based model 
can be adjusted for h > h crit to recover correctly the stationary distribution 4>{n). Thus 
we effectively replace the condition h ^> k\j [Da + Db), which requires that h is much 
larger than k\j[DA + Db) by a sharp inequality h > h crit , where h cr u is approximately 
a quarter of k\j[DA + Db)- We will show that the compartment-based model can 
be appropriately modified to correctly simulate chemical systems for any 
In particular, we also get a measure of correctness of the original compartment-based 
model. 
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3.2. Molecular-based models 

In this section, we study molecular-based models of reaction-diffusion processes, i.e. 
we simulate trajectories of individual molecules. The position [X(t),Y(t), Z(t)} of a 
diffusing molecule (Brownian motion) can be described by a system of three (uncoupled) 
stochastic differential equations (SDEs) [5] 

X(t + dt) = X(t) + V2D dW x , (15) 

Y(t + dt) = Y(t) + V2D dW y , (16) 

Z(t + dt) = Z(t) + V2D dW z , (17) 

where dW x , dW y , dW z are (uncorrelated) white noises (i.e. differentials of the Wiener 
process) and D is the diffusion constant. To simulate trajectories of the system of SDEs 
( TT5|) - (fTTl) . we choose a small time step At and use the Euler-Maruyama method to solve 
SDEs (TT5]) - (jT7j) : that is, we compute the position [X(t + At), Y(t + At), Z(t + At)] at 
time t + At from its position [X(t),Y(t), Z(t)] at time t by 



X(t + At) = X{t) + V2DAt (18) 

Y(t + At) = Y(t) + V2DAt £ y , (19) 

Z(t + At) = Z(t) + V2DAt f z , (20) 

where D is the diffusion constant and £ x , £ y , £ z are random numbers which are sampled 
from the normal distribution with zero mean and unit variance. To model a bimolecular 
reaction, it is often postulated that two molecules (which are subject to the bimolecular 
reaction) always react whenever their distance is less than a given reaction radius g 
[27] [2]. If trajectories of molecules exactly follow the system of SDEs ( fT5i) -( |T7l) . one 
can find explicit formulae linking the reaction rate constant, the diffusion constant (s) of 
reactants and the reaction radius [271 ffl E] . The reaction radius of heteroreaction (TjQ) is 

^MD^D,) (21) 
and the reaction radius of homoreaction ([5]) is k\j (8ttDa)- Thus, the illustrative example 
([!}) can be simulated as follows. We choose a small time step At. We update the 
position of every molecule by ([TBI - (120 1) where D = Da for molecules of A and D = Db 
for molecules of B. Reflecting boundary conditions are implemented on the boundary 
of the cubic computational domain [0, L] x [0, L] x [0, L]. For example, if X(t + At) 
computed by (DSD is less than 0, then X(t + At) = -X(t) - V2D At £ x . If X(t + At) 
computed by (TTg]) is greater than L, then X(t + At) = 2L-X(t)-y/2D At Similarly 
for y and ^-coordinates. Whenever the distance of a molecule of A from a molecule of 
B is less than the reaction radius g given by (pTj) . we remove the molecule of A from 
the system. We also generate a random number r uniformly distributed in (0, 1) during 
every time step. If r < k 2 L 3 At, then we generate another three random numbers r x , r y 
and r z uniformly distributed in (0, 1) and introduce a new molecule of A at the position 
[r x L, r y L, r z L). 

Considering the parameter values from Figure [3]^a), namely ki = 0.2 /im 3 sec -1 
and Da = D B = 1 /im 2 sec -1 , and using (|2~T1) . we obtain g = 8nm. On the other 
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hand, approximating the diffusing molecule as a sphere, we can estimate its radius by 
the Einstein relation [8j [3] 

e ™ = ^> < 22 > 

where k% = 1.38 x 10 -14 gmm 2 sec -2 K~ x is the Boltzmann constant, T is the absolute 
temporature, r\ is the coefficient of viscosity and D is the diffusion constant. Considering 
a solution in water (77 = 10 -3 g mm -1 sec" 1 ) at room temperature (T = 300 K), we 
obtain g m = 219.7 nm for D = Da = 1 /im 2 sec -1 . Thus the reaction radius g given 
by (|2ip is unrealistically smaller than the molecular radius g m . It is worth noting that 
this undesirable property of the model does not depend on the value of the diffusion 
constant D. For example, considering a hundred-times larger diffusion constant D, 
namely Da = Db = 100 /im 2 sec -1 , we obtain the molecular radius g m = 2.2 nm 
and the reaction radius g = 0.08 nm. Let us investigate the conditions under which 
the reaction radius g is larger than the molecular radius g m . Using (j2"2"j) . fl2Tj) and 
D = Da = Db, the inequality g > g m is equivalent to 

h > -z — • 

Considering a solution in water at room temperature, we obtain that k\ has to be at 
least of the order 10 8 M -1 sec -1 . On the other hand, typical values of k\ for interactions 
between proteins are of the order 10 6 M -1 sec -1 . Thus, unless the reaction rate constant 
k\ is very large, the model requires the reaction radius to be chosen unrealistically 
small. Notice that the values of the diffusion constant D = 1 — 100 /mi 2 sec -1 and the 
reaction rate constant ki = 0.2 /im 3 sec -1 = 3.3 x 10 6 M -1 sec -1 , which were considered 
previously, are in the range of realistic values for proteins. 

Perhaps more importantly, a small reaction radius also provides restrictions on the 
simulation time step At. We have to make sure that the average change in the distance 
between molecules during one time step, which is given by 

s = ^2(D A + D B )At, (23) 

is much less than the reaction radius g, i.e. s <C g [2J. Using Da = D B = 10 /im 2 sec -1 
and k\ = 10 6 M -1 sec -1 , we obtain that At has to be significantly less than a nanosecond. 
This limitation is even more severe for faster diffusing molecules. Note that, in the case of 
the illustrative example (jTJ), we also have to make sure that the production probability 
per one time step, k±L 3 At, is significantly less than 1. Considering the bimolecular 
reaction only, it can, in principle, be simulated with a very large time step At but the 
formula for g has to be modified accordingly. If s ^> g, then the probability that a 
given pair of molecules interacts during the time step (t, t + At) is proportional to the 
volume fraction 47r£> 3 X) / (3L 3 ) where g^ is the modified reaction radius. Comparing with 
k\ At/L 3 , we obtain 

'3k x At N 1/3 



47T 



(24) 
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This formula gives a larger (At-dependent) reaction radius, but it does not have the 
potential to provide a spatial resolution close to the size of individual molecules [2]. 
Andrews and Bray [2] designed a computational algorithm for intermediate values of At 
that satisfies s ~ g. In this case, it is not possible to derive an explicit formula relating 
g and k\ (as was done in (12~TT) for s <C g and in (T24"l) for s ^> g). Instead Andrews and 
Bray [2] provide a look-up table relating (scaled) reaction rate constant ki and reaction 
radius g. However, the reaction radius is still often smaller than molecular radius g m 
in their algorithm. In Section 14.21 we will modify molecular-based algorithms so that 
the reaction radius can be chosen as large as the molecular radius. Let us note that the 
algorithms above consider all "collisions" of reactants as reactive while in reality many 
non-reactive collisions happen before the reaction takes place. The modified algorithms 
in Section l4~2l take this point into account. 



4. Improved SSAs for reaction-diffusion modelling 

In this section, we present modified SSAs which are able to overcome the problems 
mentioned in Section [31 The make this section accessible to non-mathematicians, we 
focus only on the results. The mathematical derivation of the formulae presented and 
the justification of the modified algorithms are given in Appendices. 



4-1. Improved compartment-based model 

Let us consider the heteroreaction example (CD) modelled by the compartment-based 
react ion- diffusion model (l9])- ffTT|) . We will show that a suitable modification of 
propensity functions aijk,i(t), which were defined by ( IT2l) . leads to an algorithm that 
gives the correct <f)(n) for any h larger than or equal to the critical value h cr i t . 
Moreover, this is not possible for values of h smaller than h crit . The critical value of 
the compartment size h can be estimated as 

hcrit = Poo ~pj ; (25) 

where k\ is the rate constant of the bimolecular reaction, Da (resp. Db) is the diffusion 
constant of A (resp. B) and (3^ ~ 0.25272. If h crit computed by (1251) is significantly 
smaller than the domain size L, then the critical value of h is indeed given by ( 1251) . If 
the domain size L is comparable to (T25]) . then fT25|) provides only a good approximation 
of h cr it, with the real value being slightly higher as discussed below. If h > h cr u, we 
propose to modify the first formula in (|T2l) by 

a m At) = A it (t)B iit (t) {DA ( ^ B) _ k ; hh2 (26) 

where the parameter f3 has no physical dimension and needs to be specified. If a 
modeller does not have any information about the system, we propose to choose 
(3 = (3oo Rj 0.25272. We will show later that f3oo is indeed the correct choice of (5 if 
K = L/h is large. Formula ff26l) can be applied to any heteroreaction A + B — > where 
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K 


P 


2 


0.30208 


3 


0.33233 


4 


0.33461 


5 


0.33119 


6 


0.32660 


7 


0.32205 


8 


0.31784 



K 


P 


9 


0.31406 


10 


0.31067 


12 


0.30493 


14 


0.30027 


16 


0.29643 


18 


0.29322 


20 


0.29048 



K 


P 


25 


0.28514 


30 


0.28123 


40 


0.27587 


50 


0.27232 


60 


0.26979 


80 


0.26640 


100 


0.26420 



K 


P 


150 


0.26103 


200 


0.25930 


300 


0.25743 


400 


0.25643 


600 


0.25536 


800 


0.25479 


1000 


0.25443 



Table 1. The values of (3 for the selected values of K computed by <\27[ for the 
heteroreaction example |7]) for Bq = 1 . 



stands for an arbitrary right hand side, provided that the propensity function (126!) is 
positive. Consequently, the critical value h cr u is the one which makes the denominator 
of (126]) equal to zero. In such a case, the propensity function is infinity and the reaction 
happens immediately after the reacting molecules enter the same compartment. Thus, 
h cr it satisfies (Da + DB)h^, rit — (3h 2 crit ki = 0. If we substitute (3^ for /3, we obtain the 
approximation (1251) . 

Let us consider the illustrative heteroreaction example plotted in Figure [3](a). The 
values of the constant j3 for this model are given for different values of K in Tabled], and 
lie between (3^ rs 0.25272 and 0.34. Increasing K to infinity, the values of f3 converge to 
Poo. In Figure H](a), we compare the results computed using the original formula (fl2"|) and 
by the new formula (1261) for the heteroreaction example . We use the same parameter 
values as in Figure G^a) and K = 16. As in Figure [3](a), we observe a difference between 
4>i(n) and (fri^ri) if the original model is used. On the other hand, 4>i§(n) computed by 
the modified algorithm (solid line) is the same as <f>i(n) (grey histogram). 

In Table [Q, we observe that P weakly depends on K = L/h, so that the real value 
of h cr it (that makes the propensity function (|26l) equal to infinity) is slightly larger 
than (1251) . The dependence of P on K is caused by the boundary of the computational 
domain [0, L] x [0, L] x [0, L]. Every inner compartment can be entered from six possible 
directions, but some incoming directions are missing in the boundary compartments. 
Whenever a molecule of B is in a boundary compartment, it is less likely found 
by molecules of A. One could address this problem either by introducing different 
propensity functions for different compartments, or by using /3 in ( 1261) which is slightly 
larger than {3^. We used the latter option. In Appendix C , we show that the modified 
P is given by 

K-l 

= J— V" . (27) 

2-^ 3 i -j^L 3 — cos (iit/K) — cos (jir/K) — cos (kir/K) 

(i,j,'k)¥= (o,o,o) 
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Figure 4. (a) Heteroreaction example fl])- Stationary distribution (f>w(n) defined by 
(j 13\ computed by the original SSA (dashed line) and by the modified SSA that uses 
<\26$ instead of (ji^j) (solid line). Correct stationary distribution <fi(n) = 4>\{n) is plotted 
as the grey histogram, (b) Homoreaction example (J5)). Stationary distribution <pi6( n ) 
computed by the original SSA (dashed line) and by the modified SSA that uses <\S0$ 
instead of (solid line). Correct stationary distribution (j)(n) = 4>i(n) is plotted as 
the grey histogram. 



The results in Tabled] have been computed by (1271) . In |Appendix D[ we show that (3^ 
is given by 

i r r i 

/5oo = ^/ / ===dxdy. (28) 

^ Jo Jo a/(o — cosx — cosy) z — 1 

Evaluating this integral numerically by the Monte Carlo method, we obtain (3^ pa 
0.25272. If we model a complicated reaction-diffusion system, modelling of each 
heteroreaction will be improved by using (12"oT) . provided that all rates obtained by (|26[) 
are positive. In other words, the smallest possible h which can be simulated is given 
as the maximal h cr u for each bimolecular reaction. If h is significantly larger then h cr u 
(i.e. if h S> h crit ), then we have {Da + Db)Ii 3 3> f3kih 2 and we can approximate 
(D A + D B )k l ^ h 
(D A + D B )h 3 - f3hh 2 ~ h 3 ' 1 ; 

In particular, the propensity function Oijk,i{t) defined by ( 1261) is approximately equal 
to the original propensity function (j!2p for large values of h. On the other hand, if h 
is close to hcru, then the propensity function aiijk,i{t) given by ( |26l) is larger than the 
original propensity function (fl2l) . 

Finally, let us consider the homoreaction example (jSJ) modelled by the compartment- 
based model. In this case, we propose to modify the propensity functions (|T4l) for 
h > hcru, by a ijk>1 (t) 

Oi jk ,i(t) = Aijk(t)(A ijk (t) - 1) DA tft k p klh 2 > ( 3 °) 
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where (3 is a constant. In Figure BJb), we compare the results computed using the 
original formula (|14j) and by the new formula (130]) for the value of (3 given by Table [TJ 
We use the same parameter values as in Figure [3](b) and K = 16. We again observe 
that the modified formula (1301) gives better results than the original SSA. One can 
still observe a small error which is caused by the fact that we used the value of (3 
computed for heteroreactions. Using (1271) . we have (3 ~ 0.29643 for K = 16 (see Table 
[1]). Experimenting with the model (jSJ), we can find that (3 = 0.28 yields slightly better 
fit between 4>i(n) and 4>i6(n). Notice that f3 = 0.28 is still larger than (3^ ~ 0.25272. 
However, for the purposes of applications, it is sufficient to use either (3 = (3^ or the 
values of (3 from Table [T]for both heteroreactions and homoreactions. Using (3 = (3^, we 
discovered that the biggest contribution of the error stems from the boundary effects and 
derived Table[T]which adds a correction to to compensate for boundary behaviour. In 
a similar way, one could look for further corrections to the value of (3 for homoreactions, 
or for domains which are cuboids rather than cubes. Although such corrections are of 
interest from the mathematical point of view, they provide only a negligible improvement 
of the algorithm. Thus we will not include them in this paper. 

4-2. Improved molecular-based models 

The major assumption of molecular-based models is that molecules always react 
whenever their distance is less than the reaction radius g. The reaction radius g is related 
to the rate constant of the bimolecular reaction by a simple formula (for example, (1211) 
for the Smoluchowski model) or by a look up table for the Andrews and Bray model [2]. 
In this section, we present models that implement bimolecular reactions with the help 
of two parameters: the reaction radius g and the reaction rate A. We postulate that the 
bimolecular reaction can take place only when the distance of molecules is less than ~g. 
If this is the case, then the bimolecular reaction events happen with the rate A. We will 
call this model A — ~g model in what follows. 

To implement this idea on the computer, we need to relate the parameters A and 
~g to the rate constant of the bimolecular reaction. The advantage of the A — ~g model is 
that many different pairs of A and ~g correspond to the same bimolecular rate constant. 
In mathematical terms, the condition on A and g is one equation for two unknowns A 
and ~q. In particular, we can choose the value of ~g as desired (e.g. to be comparable to 
the molecular radius g m ) and compute the appropriate value of A. Thus A — g model has 
the potential to solve the problems of molecular-based modelling discussed in Section 

We explain the A — ~g model on the heteroreaction example ([1]). The diffusion of 
molecules A and B is simulated as in Section 13721 We choose a small time step At. The 
trajectory of every molecule is computed by (Il8l) - (j2"0l) where D = Da for molecules of 
A and D = Db for molecules of B. Let s be the average change in the relative position 
of a molecule of A and a molecule of B during one time step, given by ( 1231) . We will 
distinguish two cases of the value of the time step: (i) the time step At is chosen so 
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small that s C and (ii) the time step At is larger so that s ~ 

(%) Small time step At. To model heteroreaction A + B — > ■ ■ we introduce two 
parameters: reaction radius g and rate A. The reaction radius is expressed in units of 
length and rate A in units per time. Whenever the distance of a molecule of A from 
a molecule of B is less than the reaction radius g, then the heteroreaction takes place 
with the rate A. In Appendix E we derive the following relation between g, A and the 
rate constant k\ of the heteroreaction: 




h = 4tt(D a + D B )\-g-\l RA±£r tanh ( g \j - - - J ) . (31) 

This is one condition for two unknowns ~q and A. In particular, we can choose ~g 
comparable to the radii of reacting molecules and use (13T!) to compute the corresponding 
A. Notice that (13T1) is a simple non-linear equation which can be solved by any numerical 
method for finding roots of a real-valued function (for example, Newton's method or 
the bisection method). 

If A = oo (that is, if molecules react immediately whenever they are within the 
reaction radius), then (J3T1) simplifies to (|2T1) as desired. On the other hand, if A is small 
that A <C (Da + -Db) £? 2 , then we can use Taylor expansion in (I3~T|) to approximate 

tanh (WV(A4 + £b)) « g y/X/(D A + D B ) - l - (g y/X/(D A + D^f) \ 
Consequently, fl3T|) simplifies to k\ « 47r^ 3 A/3 which can be equivalently rewritten as 

4tt^ 3 /3' 

i.e. the reaction rate A is given as the reaction rate constant k\ divided by the volume, 
47r£> 3 /3, of the ball in which the reaction takes place. Formula ( 1321) is analogous to the 
formula for the reaction rate per compartment in the compartment-based approach for 
large compartment size h. If h is large satisfying h ^> h crit , then the reaction rate per 
compartment is given as ki/h 3 , which is the reaction rate constant k\ divided by the 
volume, h 3 , of the compartment - see (129]) . 



A « 7Z=^> (32) 



(ii) SSA for larger time steps. We introduce two parameters: reaction radius ~g and 
probability Pa- The heteroreaction A + B — > ■ ■ ■ is modelled as follows: whenever the 
distance between a molecule of A and a molecule of B (at the end of a time step) is less 
than the reaction radius ~g, then the heteroreaction takes place with probability P\; that 
is, we generate a random number r uniformly distributed in (0, 1) and the heteroreaction 
(removal /addition of molecules) is performed whenever r < P\. Notice that the previous 
algorithm (for small time step At) can be also formulated in terms of parameters g and 
P\ rather than ~g and A. Indeed, if A At 1, we have P\ ~ A At. If At is larger, then 
the relation between P\ and A is more complicated. However, from the practical point 
of view, there is no need to know the rate A: it is sufficient to formulate the algorithms 
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Figure 5. (a) Dependence of P\ on k for three different values 0/7. Dimensionless 
parameters k and 7 are given by Ijffffi . (b) Dependence of n on 7 for three different 
values of P\ . 



in terms of g and Pa- Next, we present the condition relating g and Pa with the rate 
constant k\ of heteroreaction A + B — > ■ ■ -. We define the dimensionless parameters 

s _ ^/2(D A + D B )At h At 

Q 



1 



K 



(33) 



g g" 
In applications, we first specify the time step At. We also want to specify ~q in a realistic 
parameter range. Consequently, 7 and k can be considered as given numbers in what 
follows. For example, we can choose the average change of distance between reacting 
molecules during one time step equal to the reaction radius, i.e. s = g. Then (133]) gives 
7 = 1. The key modelling question is: what is the appropriate value of the probability 
P\? In Figure [5](a), we present the dependence of P\ on k for three different values of 
7. The derivation of the equation for Pa and the numerical method which was used to 



compute this plot are given in Appendix F Below, we summarize only the equations 
that were solved and present illustrative computational results. 

To formulate the equation for Pa, it is useful to define an (auxiliary) function 
g(r) : [0, 00) — > [0, 1] as the solution of the integral equation 

pi POO 

g(r) = (1-P0 / K(r,r';j)g(r')dr'+ / K(r,r';j)g(r')dr', (34) 



satisfying g(r) — > 1 as r — ► 00, where 



K(r, r';j) 



r7v2vr 



exp 



r — r 



A2 



2 7 2 



— exp 



A2 



(r + r') 



(35) 



The function g(r) depends on dimensionless parameters P\ and 7; we make this explicit 
by writing 

9(r;Px,i) =g(r). 
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Then, the model parameters g, P\, At are related to rate constant k\ and diffusion 
constants Da, D b by 

k = P x [ 4vrr^(r;P A , 7 )dr. (36) 
Jo 

Since hi, Da and Db are known and parameters At and ~g can be specified first, 
parameters 7 and k are in applications given numbers. Thus fl36l) is one equation 



for one unknown Pa- In Appendix F we present a numerical approach for solving this 
equation, as well as the derivation of (!34|) - (l36l) . 

In Figure [5](b), we present the dependence of k on 7 for three different values of 
Pa- Note that the case Pa = 1 corresponds to the Andrews and Bray model [2]. Thus 
the solid line in Figure [S^b) has been already computed in reference [2] . However, we 
propose to use a much smaller value of Pa, which enables us to choose a larger (more 
physically meaningful) reaction radius. We see in Figure [5] that reducing P\ at 7 fixed 
corresponds to reducing k, thereby increasing the reaction radius. 

The heteroreaction example (OQ) is simulated by the A — ~g model as follows. We 
update the position of every molecule by ffT^ - fpZUl) where D = Da for molecules of A 
and D = Db for molecules of B. Reflecting boundary conditions are implemented on 
the boundary of the cubic computational domain [0, L] x [0, L] x [0, L] as in Section l3~2l 
The production of molecules of A (i.e. the second reaction in (JTJ)) is simulated as before. 
We generate a random number r uniformly distributed on [0, 1] during every time step. 
If r < k2L 3 At, then we generate another three random numbers r x , r y and r z uniformly 
distributed on [0, 1] and introduce a new molecule of A at the position (r x L,r y L,r z L). 
In particular At, has to be chosen small enough that k2L 3 At <C 1. If the separation 
between a molecule of A and a molecule of B (at the end of a time step) is less than 
the reaction radius ~g, then we generate a random number r uniformly distributed on 
[0, 1] and we remove the molecule of A from the system if r < P\. In Figure [6](a), we 
present the stationary distribution computed by this algorithm (grey histogram). The 
value of Pa, i.e P\ = 0.77%, was computed by solving (!34|) - (l36l) using the numerical 



method given in |Appendix F[ The values of parameters are given in the caption of 
Figure [6](a). Notice that the reaction radius is 40 nm and the time step At was chosen 
so that 7 = 0.5. The comparison of computational results with formula ([3]) (solid line) 
is excellent. 

Finally, let us discuss the homoreaction example (jSJ). Since two molecules of A are 
removed from the system whenever the homoreaction takes place, we have to replace 
ki by 2k\ in the above formulae. In particular, we replace k by 2k in (1361) . Moreover, 
Da + Db has to be replaced by 2D a in all formulae. Otherwise, method ( |34|) -(l36l) for 
computing P\ stays the same. In Figure E^a), we present the stationary distribution 
computed by the \-~g~ model (grey histogram). The comparison with exact formula (JTJ) 
(solid line) is again excellent. 
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Figure 6. (a) Heteroreaction example |7]). Stationary distribution computed by the \-~g~ 
model (grey histogram) and by formula ((2j) (solid line). We use ki = 0.2 /^ra 3 sec" 1 , 
&2 = 0.02 /im~ 3 sec -1 , Z?^ = Db = 10 /im 2 sec -1 , L = 2 [im and Bq = 1, 
At = 10 -5 sec, g = 40 nra and P\ = 7.7 x 10 -3 . (b) Homoreaction example (f3|) . 
Stationary distribution computed by the \-~g model (grey histogram) and by formula Q 
(solid line). We use k\ = 0.1 fim 3 sec' 1 , &2 = 0.08 fim' 3 sec' 1 , Da — 10 fim 2 sec' 1 
and L = 2 fxm, g = 40 nm and P\ = 7.7 x 10~ 3 . 



5. Discussion 



In this paper, we used the illustrative examples (OQ) and (jnj) to compare the results of 
different stochastic reaction-diffusion methods. The advantage of illustrative chemical 
models (OQ) and (jSJ) is that they have non-trivial stationary distributions given by (J3]) 
and (IT]), respectively. In principle, one could study A + B — ► B (or A + A — >• B) on its 
own to make the illustrative examples even simpler. However, the number of molecules 
of A would then decrease to zero as time progresses and the stationary distribution 
would be trivial (i.e. there would be molecules with probability 1 in the system 
after long time). The trivial stationary distribution is obtained for A + B — ► B (or 
A + A — > B) by any reaction-diffusion SSA, so we would not learn anything useful 
from the stationary behaviour. We would observe differences in modelling the transient 
behaviour of bimolecular reactions. However, the transient behaviour depends on 
the initial condition. For these reasons, we coupled bimolecular reactions with the 
production of the chemical species A to obtain the model chemical systems §T§ and (J5]) 
which have the non-trivial stationary distributions. It is worth noting that the model 
systems in this paper do not have any spatial variation of the probability distribution. 
No part of the computational domain is preferred by molecules of A or B and the 
resulting probability distribution is homogeneous in space. It is easy to generalize §B) 
and ([5]) to the spatially non-homogeneous case (for example, by considering production 
reaction — > A only in part of the computational domain p3]). Such a generalization 
is necessary for studying some other aspects of react ion- diffusion SSAs which we will 
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address in a future publication. However, our examples (JTJ and (JSJ) were complex enough 
to illustrate all results of this paper. 

We studied both on-lattice and off-lattice SSAs for reaction-diffusion processes. In 
particular, we were able to see connections between both types of models. For on-lattice 
models we found that there was a limitation on the compartment size h from below, 
i.e. h > h crit . In Section 14.11 we showed that the rate of bimolecular reaction per 
compartment must be chosen to be infinity for h = h crit . In a similar way for the 
off-lattice model, a decrease of g in the \-~g~ model, presented in Section 13.21 must be 
compensated by an increase in the rate A (probability Pa)- Again, there is a limitation 
on ~q from below, i.e. ~q must be larger than or equal to the q of the Smoluchowski model 
which is given by fl2Tj) . If ~g is sufficiently larger than this limiting value, then A is given 
by ( !32l . that is, the reaction rate constant k\ divided by the volume, 47f£ 3 /3, of the ball 
in which the reaction takes place. This is analogous to the situation h ^> h crit where 
the reaction rate per compartment is given by k\jh? (reaction rate constant k\ divided 
by the volume, h 3 , of the compartment). Thus both models give, in the limit of large ~g 
and large h, the same expression for local reaction rates: rate constant divided by the 
volume in which the reaction takes place. 

The results of this paper have been summarized in Section HI They were explained 
on illustrative computational examples, but the general formulae can be applied to 
modelling bimolecular reactions which are part of complex reaction-diffusion processes. 
We presented our results as improvements of two commonly used reaction-diffusion 
SSAs which have been previously implemented in reaction- diffusion software packages 
MesoRD [20] and Smoldyn [2]. Other molecular-based models, such as MCell [28J, 
Green's-function reaction dynamics [3.2] and velocity jump processes (TTJ [14], were 
not directly studied in this paper but some of the ideas presented in Section 14.21 can 
be applied to improve them too. From the application point of view, we focussed 
on modelling bimolecular reactions of biomolecules, e.g. proteins, but the concepts 
presented can be also applied to stochastic reaction-diffusion modelling in population 
ecology [23J or to modelling cellular dispersal [T5l 116] . In these cases, the diffusing 
objects are not macromolecules but cells or animals, and the bimolecular "reaction" is 
not a chemical reaction but local interaction between two cells or animals, for example, 
competition or predation [23] . 

Glossary 

Gillespie SSA. Stochastic simulation algorithm for simulating the time-evolution of well- 
stirred chemical systems. The results are consistent with the solution of the chemical 
master equation [TTJ EE]- 

Markov Chain. Stochastic process for which the future states of the system only depend 
on the present state and are independent of the past states. 

Modified Bessel function of the first kind. The evaluation of modified Bessel functions 
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is part of any standard mathematical software (e.g. the function besseli in Matlab). 
In general, the modified Bessel function I n (for n G N) is a solution of the ordinary 
differential equation 

z 2 r;(z) + zr n (z)-(z 2 + n 2 )i n (z) = o. 
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Appendix A. Stationary distributions, means and variances for the 
illustrative heteroreaction and homoreaction examples 

Let us consider chemical system (CQ) to be well-stirred. Let p n (t) be the probability that 
there are n molecules of A at time t in the reactor, i.e. A{t) = n. Then p n (t) evolves 
according to the chemical master equation [13], [19] 
dp n hB hBo 

— - = (n + 1) p n+ i np n + k 2 vp n -i ~ k 2 vp n (A.l) 

at v v 

where the third term on the right hand side is missing in (lA.lj) for n = 0; i.e. we use 

the convention that p_i = 0. The stationary distribution <f)(n) is defined by 

<f>{n) = limp ft (t). (A.2) 

t^oo 

Consequently, flA.lj) implies that <f>(n) satisfies the equation 

— - — - (n + 1) cj){n + 1) — -n Mn) + k 2 v 4>{n — 1) — k 2 v 4>{n) = 

v v 

where 0(— 1) = 0, which can be equivalently written as 

^ n )= (t^ + 1 -~) 0(n-l)-T^0(™-2), forn>2. (A.3) 

Thus <f)(n) is uniquely determined by the value of 0(0). We can easily verify that ([3]) 
satisfies (IA.3p . Moreover, it is the only solution of (IA.3I) that satisfies the normalization 
condition Yl'^Lo 'Pi 71 ) = 1- The stationary value of the stochastic mean M s (i.e. the 
value around which the number of molecules fluctuates) and the stationary value of the 
variance V s (i.e. the size of the stochastic fluctuations) are given by 

oo oo 

M s = ^n0(n), K = ^(n-M s ) 2 0H- (A.4) 

n=0 n=0 

Using (J3J), we obtain (jU) and V s = M s . 
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Let us consider the homoreaction example (j5J). Then the chemical master equation 
reads as follows 



dt 



— (n + 2)(n + l)p n+2 - —n(n - l)p n + k 2 vp n _ x - k 2 vp n . 



Starting with the stationary version of this equation, one can use the method of moment 
generating function [21] to show that the stationary distribution <f>(n) is given by (J7J) 
[9], [30]. The stationary values of stochastic mean M s and variance V s , which are defined 
by 0A.4p . can be also evaluated in terms of the Bessel functions; M s is given by ([S]) 
and V s = M s — M 2 + k 2 v 2 ji2k\). The classical deterministic description of the chemical 
system §5§ is given, for concentration ait) = A{t)jv, as the ODE da/dt = —2k\a 2 + k 2 . 
Multiplying by z/, we obtain the ODE 

d~A 



dt 



-2k x jvA + k 2 v 



(A.5) 



where A(t) = a{t)v is the deterministic approximation of the average number of 
molecules in the volume v with concentration ait). Notice that equation flA.51) does 
not give us the time evolution of the stochastic mean. To see that, let us consider 
the stationary value of Ait). It is given as the solution of the stationary equation 
corresponding to (1A.5I) . namely = — 2k\/vA s + k 2 i'. Hence, A s = v^Jk 2 /2ki which is 
not equal to M s given by formula (IE]). See also the discussion at the end of Section [2721 



Appendix B. Reaction-diffusion master equation 

Let N = {0, 1, 2, 3, . . .} be the set of non-negative integers. Let n e and m G W al1 . 
We denote their coordinates by three indices, namely 

{nijk I (i,3, k) E I aii} and m = {m ijk \ k) e I a u}. (B.l) 



n 



Let p(n, m, t) be the joint probability that 



and B ijk (t) 



mijk for ah 



(i,j,k) G I a ii- The reaction-diffusion master equation describes the time evolution of 
p(n, m, t). To formulate it, we define the operators J^- fc : N /a " — > N Ia " for (i,j,k) G I a ii 
and e G E ijk by 



{quvw | {U,V, W) G I a ll} 



where 



^UVW 5 



+ 1, 
" 1, 



for (w, v, w) 
for (w, v, w) 
otherwise. 



k) +e; 



We also define 



'ijk 



{5 UVW | (u, v, w) G I a ii} where S 



1, for (w, v, w) 
0, otherwise. 



k); 
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Then the reaction-diffusion master equation, i.e. the chemical master equation which 
corresponds to the system of "chemical reactions" fl9l)- f|TTl) . can be written as follows 

P21EEE] 

dp(n, m) k\ 



dt 

{i,j,k)£l a u 



r| Y {( n ijk + l)^ufcP(n + $ijk, m) - n ijk m ijk p(n, m)| 

(i,j,k)eI aU 

+ k 2 h 3 Y [p(n-S ijk ,m) -p(n,m)| (B.2) 

(ij,/c)G/aii 

+ "^T E E {( n ^ + 1 M J ^( n ), m ) -n^p(n,m)J 
+ "^f E E {( m ^ + ^ p ( n ' J 5fc( m )) ~ m iifcP( n > m )}- 

(i,j,k)el a u eSEijfe 

The first two terms on the right hand side correspond to chemical reactions (Q, the 
third term to diffusion jumps ffTUl) and the last term to diffusion jumps ffTTl) . 

Appendix C. Derivation of formulae ( I26|) and ( 1271) 

We will first study the case Db = and £?o = 1. This means that there is only one 
molecule of B in the system and it does not diffuse. In particular, reactions (II ip are not 
included in the model. Let the molecule of B be in the compartment b — (bi, &2, 63) G ^ozz- 
Let n e N /a " with the coordinates defined by (IB. II) . Let p(n,t) be the joint probability 
that Aijkit) = n^k for all (i, j, k) e J a ^. Since the position of the molecule of B does not 
evolve, the reaction-diffusion master equation (IB. 21) simplifies to the following equation 
for p(n, i) 

= ^{(%+ !)p( n + <%) - %K n )} + ^ 3 5^ {p(n - S ijk ) - p(n)} 

(i,j,k)el a ii 

+ E E {K* + i)pW n ))-%K")}. (ci) 

(i,j,k)ei aU eeE^ fc 

We want to change the reaction rate ki/h 3 (of the bimolecular reaction per one 
compartment) to a reaction rate A in order to decrease the error between stationary 
distributions 4>k and 0i. The stationary version of (IC.lj) with ki/h 3 replaced by A is 



^{( n b+ 1 )Ps(n + (5 F ) -7ijp a (n)| + fc 2 /i 3 |p s (n - <5y fc ) -p s (n)| 

(i,j,k)el a u 

+ ^ E E {(^ + l)Ps(^ fc (n))-n ufcPs (n)} = (C.2) 



(i,i,fe)e/ « eeE ijfe 

where 



p s (n) = lim p(n, t). (C.3) 

t— >oo 
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Let us denote the average number of molecules at the lattice site k) as 



oo oo 



M ijk (t) = y^njjfcp s (n) = ^ ••• ni, fc p s (n) 

n "000=0 nooi=0 "ifif/f=0 



Multiplying (1C2|) by riy-jt and summing over n, we obtain 

h 2 



k 2 h 3 + 2A. (Mi, k+e - M ijk ) = 0, for k) ± b, (C.4) 
k ^ + TT E ( M M- - = AM,. (C.5) 



eeE F 



Let us define tensors ^ j k G ]R^x^x^ ; for i' = 0, 1, 2, . . . , K - 1; / = 0, 1, 2, . . . , K - 1 
and k' = 0, 1,2,. . .,K - 1, by 

W 1 f i'ji -l/2)n \ f f(j - l/2)n \ fk'(k - l/2)vr 

V'i/fc = cos x cos x cos 



\/8, if i', j', k! are nonzero; 

V4, if exactly one of i', j', /c' is zero; 

\/2, if exactly two of z', j', fc' are zero; 

1, for i! = f = k! = 0. 



x < 



(C.6) 



We have 

K 



E cos — cos i/ = for 1 > °' 



K I V K / 2 

where is the Kronecker delta. Consequently, ip l 'i' k ' , for i', j', fc' = 0, 1, . . . , K — 1, 
satisfy the ort honor mality condition: 

K 

E ^i/^ ''4jk k " = 5 i'i" S fj" S k'k"- (C7) 
i,j,fe=l 

Let us express My k in the basis il) % '^' k ' : 

K-l 

M l3k = Y, M i/j/k ^% k '. (C.8) 

i',j',k'=0 

Then (IU4]l -(lC31) read as follows 

D A 11 



k ^ + J^ E Mi>fk> £ =0. for (<, J,*) ^5, 

i',j',k'=0 eeE iifc 
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Using (10.61) . we get 

n ~ 

k2Ki + li E Mi'fk> c i>j ' k ' ij}\jl k = 0, for(i,j,fc)^6 l (C.9) 



where 



i',j',k'=0 

k2h 3 + ^ M iW /^>f fc ' = AM F) (CIO) 

i',j',fc'=0 

c .^ 2 (cos(|) + cos(f) + cos(^)-3). (C.ll) 

Multiplying flC.9p - flC.10p by fc , summing the resulting equations and using the 
orthonormality condition (1C.7I) . we obtain 

k 2 v = AA%, (C.12) 



^^"*"M 4W = AM^f for (i",j",k") + (0,0,0), 



h 2 

where v = L 3 = h 3 K 3 is the volume of the reactor. We drop the double primes on 
indices i, j and k to simplify the notation and obtain 



D A c^ k 



Mijk = - 'J , for (i,j,k) ? (0,0,0). 
Using fiOT2|) . we get 



k v ib^ k 

M ijk = \ 2 , for (i,j,k)^ (0,0,0). (C.13) 



Using (JUp and flCTT3l) . we have 



ij,fe=0 

where 

-ft"— 1 C„i,ijk\2 



A 



(ib- , 

/%=- E t ( c - 15 ) 



ijjjfc = 
(i,j,fc) ^ (0,0,0) 

Substituting fRXl2l) into fjCHIj) . we have 

^^-"'M™-^. (C.16) 
The average number of molecules of A in the reactor is given by 

K 

A s = M Hk- 

i,j,k=l 

Using (1C.8p . we get 

^ = E E ^'i'*'V>g fc ' = E ^' = Mooo^ 3/2 . (C.17) 

i,j,k=l i' ,j' ,k'=0 i',j',k'=0 i,j,k=l 
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We would like to choose A so that A s = M s where M s is given by (j3J, i.e. 
h " 2 M 000 K^. 



h 

Substituting for M 000 into (1C.16H and using v = h 3 K 3 , we get 

, * kov „ kov h 2 



A k\ Da 

which implies 

A _ D A h 

D A h 3 - fchh 2 ' 

This choice of A gives the average number of molecules of A equal to M s , provided that 
Db — 0, Bq — 1 and the molecule of B is in the compartment b = (&i, 62, 63) £ laii- 

Now let Db 7^ and B — 1. If we want to model the bimolecular reaction (JTj), it is 
important to know the distances of molecules of A from the molecule of B. The distances 
diffuse with the diffusion constant Da + D B - Thus we can equivalently model their time 
evolution by considering that the molecule of B does not diffuse and molecules of A 
diffuse with the diffusion constant Da + D B - Then the previous calculation (equation 
([CT6]) ) implies that 

K^M^ l ^ + lh J^ (C.18) 

A L) a + i->B 

where notation Mq 00 = M 000 highlights the fact that M 000 depends on the position b of 
the molecule of B. Let be the probability that the molecule of B is in the compartment 
b = (bi,b 2 ,b 3 ) G I aii- We have p^ = K~ 3 . The average number of molecules of A in the 
reactor is given by (compare with (1C.17I) ) 

b b 
Multyplying (jC18l) by p^ = K~ 3 and summing over 6, we obtain 

b b 

We would like to choose A so that A s = M s where M s is given by (J1J. Substituting (jlj) 
into (1C.19j) . we get 

k-yv 2 kov „ kovh 2 



where 



k x K 3 A ^ Da + Db 
Using v = h 3 K 3 , we obtain 

{D A + D B )h 



A 



(Da + Ds)h 3 — (3k\h 2 
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Thus we have derived (1261) . Using (lC.20p . (1C.15P and orthonormality condition (lC.7p . 
we obtain 

1 1 K ~ x (ib^) 2 1 K ' x 1 

8=— V/3-=-— V V b - V — 

~b i,j,k = Jy i,j,k = 

(i,j,k) + (0,0,0) (i,j,k) + (0,0,0) 

Substituting flCTTTTh for c ij ' fc , we obtain 1127) 1 . 
Appendix D. Derivation of formula ( 128ft 



Formula (127)) is the Riemann sum of the definite integral 

1 r 7 ^ r 7 ^ z* 71 " 1 

p oo = -^ / / dxdydz, (D.l) 

2 ( 7r ) Jo Jo Jo 3 - cosx - cosy - cosz 



i.e., passing .K" — > oo in ( l27j) . we obtain ( ID. lh . Integrating over z, we get ( 1281) . 

Formula (128]) can be also derived directly without the help of (127)) . Such a derivation 



uses a similar reasoning as the derivation of (|3T|) in Appendix E, i.e. it establishes a link 
between molecular-based models and the compartment-based modelling. We consider 
the infinite three-dimensional lattice 

k)h for i G Z, j G Z, fc G Z (D.2) 

where /iGl. To model bimolecular reactions by the compartment-based approach, we 
need to know whether the molecules are in the same compartment or not. In particular, 
it is sufficient to track the relative distance of molecules rather than their absolute 
positions. Postulating that the molecule of B is always at the origin (compartment 
(0,0,0)) and letting the molecule of A diffuse with the diffusion constant (Da + Db), 
we obtain the stochastic model which gives the same distribution of relative distances 
of molecules as the original stochastic model. Thus we will study the following auxiliary 
stochastic process. We consider that the particles jump to neighbouring lattice sites 
with the rate (Da + Ds)/h 2 and are removed at the origin with the rate A. The 
react ion- diffusion master equation can be written for this model as follows (using the 
same notation as in (]B.2j) ) 

dp(n) D A + D B 



Y Y + 1 )p(« / y*( n )) _ n vkP(n)} 



dt h 2 

{i,j,k)eZ 3 eGE 

+ A|(n 00 o + l)p(n+ S 000 ) - n 00 op(n)|. (D.3) 

We are interested in the stationary behaviour of a system of (infinitely) many molecules 
of A, subject to the condition that the average number of molecules per compartment is 
kept constant (equal to M^) far from the origin, i.e. in the limit \/i 2 + j 2 + k 2 — > oo. 
The stationary version of (1D.3|) reads as follows 
D A + D B 



h 2 

(ij',fc)£Z 3 eeE 



Y Y + Ps( J ijk( n )) - n ijk p s (n)j 
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= -AjOooo + l)Psfa + ^ooo) - ^oooPs(n)| (D.4) 

where p s (n) is defined as in ( 1C.3I) . Let us denote the average number of molecules at 
the lattice site k) as 

M ijk (t) = 2jny fc p s (n). 

n 

Multiplying (1D.4|) by and summing over n, we obtain 
D A + D B 



h 2 

eeE 



J2( M m+ e - M ijk ) = for (i, j, k) ± (0, 0, 0), 



eeE 

Let us define = — M^. Then we have 

^ Mijfc+e = GfMjk, for (z,j, fc) ^ (0, 0, 0), 

eeE 

A ft 2 

V" /i e = 6yU 00 + 7^ — pr~ (/"OOO + Mx.) • 

Multiplying by exp im exp 1 ^ exp lzk , where i = \J— 1, i 6 K, 1/ 6 1, z 6 K, and summing 
over z, j and fc, we obtain 

6 % yz = % yz (exp lx + exp~ lx + exp iy + exp~ iy + exp lz + exp~ lz ) 

-\h 2 (/iooo + Moo) /{D A + D B ) (D.5) 

where Ji xyz is the Fourier transform 



oo oo oo 



%yz=Yl Yl ^P ixi exp iyj exp izk fl ijk . 

i=—oo j=—oo k=—OQ 

Simplifying (ID .5j) . we obtain 

A h 2 (^ooo + Moo) 1 



2(Da + D b ) cosx + cosy + cos z — 3 
Thus 



A h 2 (^ooo + Mqq) mR v 
/^ooo = -Poo p, — (D.6) 

L>a + L>b 



where /?oo is the constant given by 



I p2w j-2-K i-2-k | 



/5oo = — — / / / dxdy&z. (D.7) 

2(27r)' :i J J J 3 — cosx — cosy — cos z 

Using /i oo = M oo — Mx., equation flD.61) can be rewritten as 

The term AM oo gives the rate of removal of molecules of A at the origin. The rate 
of change of the concentration a of molecules of A, which is subject to heteroreaction 
([[]), can be also described by the deterministic ODE da/dt = —k\ab where b is the 
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concentration of molecules of B. This ODE can be equivalently rewritten in terms of 
the average numbers of molecules of A and B per lattice site, i.e. in terms of A = ah 3 
and B = bh 3 , as dA/dt = —k\/h 3 AB. Using B = 1 and A = Moo, the rate of removal 
of molecules of A is given by kxjh^M^. Comparing with (ID. 81) . we obtain 



h A (Da + D B ) 

h 3 00 (D A + D B ) + \{3 00 hi 



Solving for A, we obtain 



A 



(D A + D B )h 



(Da + DbW-P^W 

Thus we have derived ( |26i) with (3 = (3^. The constant (5^ is given by ( ID. 71) . Using 
periodicity of the cosine function, we obtain (ID.lj) . Integrating over z, we derive ( ]28l) . 



Appendix E. Derivation of formula (1311) 

In order to derive (l3"Ti) , we consider the diffusion to the ball of radius ~g which removes 
molecules of A with the rate A. Let the centre of the ball be at the origin. Let c(r) be 
the equilibrium concentration of molecules of A at distance r from the origin, which is 
a continuous function with continuous derivative satisfying the equations 



d 2 c 2 dc 

dr 2 r dr 

d 2 c 2 dc 
+ — r- 



Ac 



0. 
0. 



for r > g, 
for r <~g. 



dr 2 r dr Da + Db 
The general solution of these second-order ODEs can be written in the following form 

a 2 



c(r) = ai + 



c(r) 



o 3 



exp 



r 



A 



Da + D 



B 



CZ4 

H exp 

r 



-r 



A 



A4 + £ 



for r > g, 
for r <~g, 



where ai, 02, 03 and 04 are real constants. We impose the boundary condition at infinity 
lim c(r) = Coo. 

This implies a\ = Coo- Since c is continuous at the origin, we deduce 04 = — a 3 . Thus 
we have 



c(r) 



cir 



02 

CoO T , 

r 



3 sinh I r , 
/• \ \J D A + D B 



A 



for r > g, 
for r < 7j. 



To determine the constants 02 and 03, we use the continuity of c and its derivative at 
r = ~g. We obtain 

a 2 
fl.3 



/ (D A + J D B )/A tanh (g y/\/(D A + D B j) - 7i} 

Coo V(Da + D b )/\ (2 cosh (t;v / V(^a + j Db) 
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The flux through the unit area of the boundary can be computed as 

(D A + D B )a 2 



The area of the sphere is 4ng 2 . Thus the total flux through the sphere boundary is 
—Anx[D A + D B )a 2 . Substituting for a 2 , we get 

4ir(D A + D B ) (g - ^(D A + D B )/\ tanh (g y/\/(D A + D B ty Coo . 

This quantity is equal to the rate constant of bimolecular reaction k\ multiplied by 
the concentration of the chemical far from the reacting molecule Cqo. Dividing by Cqo, 
we derive (I3ip . Let us note that we used diffusion to the ball to derive (1311) . This 
approximation can be justified using the more general evolution equation for the many 
particle distribution function [7j. 

Appendix F. Derivation of (134R — (136|) and a numerical method for solving it 

Let Ci(r) be the concentration of molecules of A at distance r from the origin. Assuming 
that molecules of A only diffuse, their concentration at point r after the time interval 
At is given as 

/•oo 

K(r, r'; 7) Cj(r') dr' (F-l) 







where K(r,r'] r y) is given by (!35j) . Let us assume that the particles are removed, in the 
circle of radius g and centered at origin, with probability P\, and then diffuse for time 
At. Then ( 1F.1I) is modified to 



a l+1 {r) = {1 - P x ) J A(r,r';7) Ci (r')dr' + y K(r, r'; 7) air') dr'. 

Equation ( 134|) is an equation for the fixed point of this iterative scheme. The function 
g(r) is the generalization of the radial distribution function (RDF) for bimolecular 
reaction at steady state [2] for arbitrary Pa G [0, 1]. Note that the RDF in [2] was 
only computed for P\ — 1. The rate of removal of particles (at steady state) during one 
time step is given by the right hand side of ( 136|) . Comparing with k, we obtain ( T36l ). 

To solve (1341) . we will use the condition g(r) — >• 1 as r — > 00. Choosing S large, 
we can approximate g(r) = 1 for r > S. Let N\ and N 2 be positive integers. We 
consider the mesh rj = j/N\, for j — 1, 2, . . . , Ni and r 3 - = 1 + (S — l)(j — Ni)/N 2 , for 
j = N x + 1, . . . , Nt + N 2 . We discretize (El as 



Y_ p N S - I Nl+m r°° 

gin) = N X K ( ri > r i> ^9(rj) + -ji— K ( ri ' r i ; ^9(rj) + / K(r it r'; 7) dr'. 

This is a linear system for g(ri), i — 1, 2, . . . , N± + N 2 , which can be solved, for example, 
by Gaussian elimination. Let us note that the right hand side of this system can be 
evaluated using the error function erf as 



K{r h r ;7)dr = h 1 - - erf 



1V2 



- erf 
2 



S + r t 
7V2 



Stochastic reaction- diffusion processes 



32 



Substituting gfa), i = 1, 2, . . . , Ni + N 2 , into (13T)j) . we compute n. Repeating this 
computation for different values of 7 and P\, we obtain the results presented in Figure 

El 
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